home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / FORTRAN1.LZH / INTRPL.FOR < prev    next >
Text File  |  1988-02-08  |  4KB  |  166 lines
  1.       SUBROUTINE INTRPL ( L, X, Y, N, U, V, IERR )
  2. C
  3. C  INTRPL INTERPOLATES A FUNCTION Y(X) FROM A SET OF VALUES MONOTONICALY
  4. C  INCREASING IN X.
  5. C
  6. C  INPUTS :
  7. C    L = NUMBER OF DATA POINTS (>= 2)
  8. C    X = ARRAY OF X VALUES IN INCREASING ORDER
  9. C    Y = ARRAY OF Y VALUES
  10. C    N = NUMBER OF POINTS TO BE INTERPOLATED
  11. C    U = ARRAY OF X LOCATIONS FOR OUTPUT (N VALUES)
  12. C
  13. C  OUTPUTS :
  14. C    V = ARRAY OF INTERPOLATED Y VALUES
  15. C    IERR = ERROR INDICAT0R...
  16. C             0 = NO ERROR
  17. C             1 = L TOO SMALL
  18. C             2 = N LESS OR EQUAL TO 0
  19. C             3 = X VALUES NOT MONOTONICALY INCREASING
  20. C
  21.       DIMENSION X(L), Y(L), U(N), V(N)
  22.       REAL M1,M2,M3,M4,M5
  23.       EQUIVALENCE (X2,A1,M1),(X5,A5,M5),(Q1,T3)
  24.       EQUIVALENCE (J,SW,SA),(Y2,W2,W4,Q2),(Y5,W3,Q3)
  25. C
  26.       IERR = 0
  27.       LO  = L
  28.       LM1 = LO - 1
  29.       LM2 = LM1 - 1
  30.       LP1 = LO + 1
  31.       NO  = N
  32.       IF ( LM2 .LT. 0 ) THEN
  33.          IERR = 1
  34.          RETURN
  35.       ENDIF
  36.       IF ( NO .LE. 0 ) THEN
  37.          IERR = 2
  38.          RETURN
  39.       ENDIF
  40.       DO 11 I = 2, LO
  41.          IF(X(I-1) .GE. X(I)) THEN
  42.             IERR = 3
  43.             RETURN
  44.          ENDIF
  45. 11       CONTINUE
  46.       IPV = 0
  47.       DO 80 K = 1, NO
  48.          DX = U(K)
  49.          IF (LM2 .EQ. 0) THEN
  50.             I = 2
  51.          ELSE IF (DX .GE. X(LO)) THEN
  52.             I = LP1
  53.          ELSE IF (DX .LT. X(1)) THEN
  54.             I = 1
  55.          ELSE
  56.             IMN = 2
  57.             IMX = LO
  58. 21          I = (IMN + IMX)/2
  59.             IF (DX .LT. X(I)) THEN
  60.                IMX = I
  61.             ELSE
  62.                IMN = I + 1
  63.             ENDIF
  64.             IF (IMX .GT. IMN) GO TO 21
  65.             I = IMX
  66.          ENDIF
  67.          IF (I .EQ. IPV) GO TO 70
  68.          IPV = I
  69.          J = I
  70.          IF (J .EQ. 1)J=2
  71.          IF (J .EQ. LP1)J=LO
  72.          X3 = X(J-1)
  73.          Y3 = Y(J-1)
  74.          X4 = X(J)
  75.          Y4 = Y(J)
  76.          A3 = X4 - X3
  77.          M3 = (Y4 - Y3)/A3
  78.          IF (LM2 .EQ. 0) THEN
  79.             M2 = M3
  80.             M4 = M3
  81.             GO TO 45
  82.          ENDIF
  83.          IF (J .NE. 2) THEN
  84.             X2 = X(J-2)
  85.             Y2 = Y(J-2)
  86.             A2 = X3 - X2
  87.             M2 = (Y3 - Y2)/A2
  88.             IF (J .EQ. LO) THEN
  89.                M4 = 2*M3 - M2
  90.                GO TO 45
  91.             ENDIF
  92.          ENDIF
  93.          X5 = X(J+1)
  94.          Y5 = Y(J+1)
  95.          A4 = X5-X4
  96.          M4 = (Y5 - Y4)/A4
  97.          IF (J .EQ. 2) M2 = 2*M3-M4
  98. 45       IF (J .LE. 3) THEN
  99.             M1 = 2*M2 - M3
  100.          ELSE
  101.             A1 = X2 - X(J-3)
  102.             M1 = (Y2-Y(J-3))/A1
  103.          ENDIF
  104.          IF (J .GE. LM1) THEN
  105.             M5 = 2*M4 - M3
  106.          ELSE
  107.             A5 = X(J+2) - X5
  108.             M5 = (Y(J+2) - Y5)/A5
  109.          ENDIF
  110. C
  111. C --- NUMERICAL DIFFERENTIATION
  112. C
  113.          IF (I .EQ. LP1) GO TO 52
  114.          W2 = ABS(M4-M3)
  115.          W3 = ABS(M2-M1)
  116.          SW = W2 + W3
  117.          IF (SW .EQ. 0) THEN
  118.             W2 = .5
  119.             W3 = .5
  120.             SW = 1.0
  121.          ENDIF
  122.          T3 = (W2*M2 + W3*M3)/SW
  123.          IF (I .EQ. 1) GO TO 54
  124. 52       W3 = ABS(M5-M4)
  125.          W4 = ABS(M3-M2)
  126.          SW = W3 + W4
  127.          IF (SW .EQ. 0.) THEN
  128.             W3 = .5
  129.             W4 = .5
  130.             SW = 1.0
  131.          ENDIF
  132.          T4 = (W3*M3 + W4*M4)/SW
  133.          IF (I .EQ. LP1) THEN
  134.             T3 = T4
  135.             SA = A2 + A3
  136.             T4 = .5*(M4 + M5 + A2*(A3-A2)*(M2-M3)/SA**2)
  137.             X3 = X4
  138.             Y3 = Y4
  139.             A3 = A2
  140.             M3 = M4
  141.          ENDIF
  142.          GO TO 60
  143. 54       T4 = T3
  144.          SA = A3 + A4
  145.          T3 = .5*(M1+M2+A4*(A4-A3)*(M3-M4)/SA**2)
  146.          X3 = X3 - A4
  147.          Y3 = Y3 - M2*A4
  148.          A3 = A4
  149.          M3 = M2
  150. C
  151. C --- DETERMINATION OF COEFFICIENTS
  152. C
  153. 60       Q2 = (2.0*(M3-T3)+M3-T4)/A3
  154.          Q3 = (T3+T4-2.*M3)/A3**2
  155. C
  156. C --- COMPUTE POLYNOMIAL
  157. C
  158. 70       DX = DX - X3
  159.          V(K) = Y3+DX*(Q1+DX*(Q2+DX*Q3))
  160. 80       CONTINUE
  161.       RETURN
  162.       END
  163. C
  164. C---END INTRPL
  165. C
  166.